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Abstract 

We discuss the symbolic dynamics of biochemical networks with separate 
timescales. We show that symbolic dynamics of monomolecular reaction net¬ 
works with separated rate constants can be described by deterministic, acyclic 
automata with a number of states that is inferior to the number of biochemical 
species. For nonlinear pathways, we propose a general approach to approximate 
their dynamics by finite state machines working on the metastable states of 
the network (long life states where the system has slow dynamics). For net¬ 
works with polynomial rate functions we propose to compute metastable states 
as solutions of the tropical equilibration problem. Tropical equilibrations are 
defined by the equality of at least two dominant monomials of opposite signs 
in the differential equations of each dynamic variable. In algebraic geometry, 
tropical equilibrations are tantamount to tropical prevarieties, that are finite 
intersections of tropical hypersurfaces. 


1 Introduction 

Networks of biochemical reactions are used in computational biology as mod¬ 
els of signaling, metabolism, and gene regulation. For various applications it 
is important to understand how the dynamics of these models depend on in¬ 
ternal parameters and environment variables. Traditionally, the dynamics of 
biochemical networks is studied in the framework of chemical kinetics that can 
be either deterministic (ordinary differential equations) or stochastic (continu¬ 
ous time Markov processes). Within this framework, problems such as causality, 
reachability, temporal logics, are hard to solve and even to formalize. Concur¬ 
rency models such as Petri nets and process algebra conveniently formalize these 
questions that remain nevertheless difhcult. The main source of difficulty is the 
extensiveness of the set of trajectories that have to be analysed. Discretisation 
of the phase space does not solve the problem, because in multi-valued networks 
with m levels (Boolean networks correspond to to = 2) the number of the states 
is to" and grows exponentially with the number of variables n. An interesting 
alternative to these approaches is symbolic dynamics which means replacing 
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the trajectories of the smooth system with a sequence of symbols. In certain 
cases, this could lead to relatively simple descriptions. According to the famous 
conjecture of Jacob Palis m, smooth dynamical systems on compact spaces 
should have a finite number of attractors whose basins cover the entire ambient 
space. Compactness of ambient space is satisfied by networks of biochemical re¬ 
actions because of conservation, or dissipativity. For high dimensional systems 
with multiple separated timescales it it reasonable to consider the following 
property: trajectories within basins of attraction consists in a succession of fast 
transitions between relatively slow regions. The slow regions, generally called 
metastable states, can be of several types such as attractive invariant manifolds, 
Milnor attractors or saddles. Because of compactness of the ambient space and 
smoothness of the vector fields defining the dynamics, there should be a finite 
number of such metastable states. This phenomenon, called itinerancy received 
particular attention in neuroscience |18j . We believe that similar phenomena oc¬ 
cur in molecular regulatory networks. A simple example is the set of bifurcations 
of metastable states guiding the orderly progression of the cell cycle. In this 
paper we use tropical geometry methods to detect the presence of metastable 
states and describe the symbolic dynamics as a finite state automaton. The 
structure of the paper is the following. In the second section we compute the 
symbolic dynamics of monomolecular networks with totally separated constants. 
To this aim we rely on previous results miiiisi- In the third section we intro¬ 
duce tropical equilibrations of nonlinear networks. Tropical equilibrations are 
good candidates for metastable states. More precisely, we use minimal branches 
of tropical equilibrations as proxys for metastable states. In the forth section 
we propose an algorithm to learn finite state automata defined on these states. 

2 Monomolecular networks with totally sepa¬ 
rated constants 

Monomolecular reaction networks are the simplest reactions networks. The 
structure of these networks is completely defined by a digraph G = {V,A), in 
which vertices iGl/,l<i<n correspond to chemical species Ai, edges (z, j) S 
A correspond to reactions Ai —>■ Aj with kinetic constants kji > 0. For each 
vertex, Ai, a positive real variable Ci (concentration) is defined. The chemical 
kinetic dynamics is described by a system of linear differential equations 

( 1 ) 

3 3 

where kji > 0 are kinetic coefficients. In matrix form one has : c = Kc. The 
solutions of 0 can be expressed in terms of left and right eigenvectors of the 
kinetic matrix K: 

n—1 

c(t) = c(0)) + c(0)) exp(Afet), (2) 

k=l 

where A, are right and left eigenvectors of K, KA = XkA, = \kl^. 

The system 0 has a conservation law ^(ci-|-C 2 -|-.. ACn) = 0, and therefore 
there is a zero eigenvalue Aq = 0, = (1, 1,..., 1), (Z°, c(0)) = ci(0) -I- C2(0) -I- 

... -I- c„(0). We say that the network constants are totally separated if for all 
(z,j) ^ (z^/) one of the relations kji A. kj'i', or kji A kj>i> is satisfied. 
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It was shown in [U [T^J [13] that the eigenvalues and the eigenvectors of an 
arbitrary monomolecular reaction networks with totally separated constants can 
be approximated with good accuracy by the eigenvalues of and the eigenvectors 
of a reduced monomolecular networks whose reaction digraph is acyclic (has no 
cycles), and deterministic (has no nodes from which leave more than one edge). 
Let us denote by Gr = (Vr^Ar) the reduced digraph, and by Ki the kinetic 
constant of the unique reaction that leaves a node i € Vr- The algorithm to 
obtain G from Gr can be found in [Ida uni and will not be repeated here. 
Because Gr is deterministic it defines a flow (discrete dynamical system) on 
the graph: $(i) = j, where j is the unique node following i on the digraph. 
Reciprocally, we define Pred(i) = as the set of predecessors of the node 

i in the digraph Gr, namely Pred(i) = {j S Vr\{j,i) G Ar}- 

We say that a node is a sink if it has no successors on the graph. For the 
sake of simplicity, we suppose that there is only one sink. For each one of the 
remaining n — 1 nodes there is one reaction leaving from it. For a network with 
totally separated constants we have 

Ki Kj, or Hi 3> Kj for all i, j G [l,n — 1], i ^ j (3) 

For totally separated constants the following lemma is useful 


Lemma 2.1 If ([^ is satisfied then, at lowest order, we have 


Kj 

-Kk + Kj 


1; ifi= j o.nd Kk < Ki 

— 1, if i = k and Kj < Ki 

0, if Ki < V[\iTL{Kk, Kj) 
±oo, else 


(4) 


The dynamics of the reduced model is given by 

(5) 

jGPred(2) 

where Pred(i) is the set of predecessors of the node i in the digraph Gr, namely 
Pred(z) = {j G Vr\{j,i) G Ar}. 

As shown in ^ the eigenvectors of the approximated kinetic matrix satisfy 



= (A -1- Kfiri 

(6) 

jGPred(z) 




= (A -f Ki)li, 

(7) 


where A is the eigenvalue, Vi, h, 1 < i < n are the components of the right and 
left eigenvectors, respectively. 

Eqs. and Q imply that the right and left eigenvectors can be computed 
by recurrence on the graph, in the direct direction and in the reverse direction, 
respectively. In order to have non-zero eigenvectors, A = —Ki for some i not a 
sink, therefore the (non-zero) eigenvalues are \k = —Kk, I < fc < n — I. Taking 
into account the separation conditions ^ we get the following 


Proposition 2.2 Let us consider that Kk = 0 when k is a sink in the graph Gr- 
Then, the eigenvalues of the kinetic matrix with totally separated constants are 
Afe = —Kk, with Xk = 0 when k is a sink. The corresponding left eigenvectors 
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are 


jk _ } 1 ) if = k for some m > 0 and > Kk for all I = 0,... ,m — 1 


0 , otherwise 


and the right eigenvectors are 
1 , if j = k 

k _ \ if j = for some m > 0 and < Hit><-{k)^ 

I for all I = 1,... ,m — 1 


( 8 ) 


0 , 


otherwise. 


The full proof of the Proposition |2 . 2| can be found in the appendices. 

Let us now discuss the symbolic dynamics of the system. For each eigenvalue 
Afe = — Kfe, Kfc > 0 we associate a transition time tk = . Without loss of 

generality we can consider that ti <C ^2 ^ ^ tn-i- Any trajectory of the 

system is given by (i). At the time one exponential term exp{\kt) will vanish 
and the result will be a transition c —>■ c — r^[l^, c( 0 )), provided that {l^, c( 0 )) 

0. In other works, a trajectory can be described as a discrete sequence of 
states c(0),c(0) — r^(Z^,c(0)),.... Let us consider the following normalization 
ci(0) + 02 ( 0 ) + ... + c„(0) = 1. Then a is the probability of presence in the 
node i of a particle moving through the reaction network. For monomolecular 
networks, particles are independent, therefore this simple picture is enough for 
understanding the dynamics. Let the index iq define the initial state of the 
system Ci(,(0) = 1, Cj(0) = 0 for j 7 ^ ip. iq represents the initial position of 
the particle. According to the Prop. 2.2 (Z^,c(0)) = I^q = 1 if the step ku is 


downstream of *o in the graph and if all steps from zg to k are faster than Kk ■ 
In this case the jump at tj, is —r^. A jump —r^ has two components different 
from zero, — r* = —1 and —rj = 1 , where j is the first node downstream of 
k from which starts a step slower than Kk- Thus, the jump —corresponds 
to displacing the particle from k to j. The set of right eigenvectors defines a 
symbolic flow on the reaction digraph. A particle starting in zg first jumps in 
Zi where Zi is the first node such that < Kig, then continues to Z 2 where 
Z 2 is the first node such that Ki^ < and so one and so forth until it gets 
to the sink. Some nodes have negligible sojourn time, namely nodes such that 
Ki > Kj for all j € Pred(z). This proves the main result of the section. 

By transition graph of a finite state machine we mean the digraph Grs = 
(Vs, As), where Vs is the set of states of the machine and (i,j) € As if there are 
transitions from the state i to the state j. We have the following theorem: 


Theorem 2.3 The symbolic dynamics of a monomolecular network with totally 
separated constants can be described by a deterministic acyclic finite state ma¬ 
chine. The transition graph Grs = (14) As) of this machine can he obtained 
from the graph Gr = (Vr,Ar) in the following way: 14 = W \ {* G Vrl^i > 
Kj for all j G Pred(z)}, As = {{i,j)\i,j G 14 and there are io = i,ii,... ,im = 
j, such that ii S 14 \ 14) fori = 1,..., to — 1, and (z/)Z;+i) G Ar fori = 
0,..., TO — 1}. 

Remark 2.4 An example is detailed in Figure^ 
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Figure 1: Symbolic dynamics of a monomolecular network with total separa¬ 
tion. The integers 7 ^ labelling the reactions represent the orders of the kinetic 
constants, smaller orders meaning faster reactions. The model was reduced 
using the recipe described in HIS] (see appendices), a) full model; b-c) re¬ 
duced model with active transitions and corresponding eigenvectors. During a 
transition the network behaves like a single step : the concentrations of some 
species (white) are practically constant, some species (yellow) are rapid, low 
concentration, intermediates, one species (red) is gradually consumed and an¬ 
other (pink) is gradually produced. The net result is the displacement of a 
particle one or several steps downstream; d) The transition graph of the finite 
state machine representing the symbolic dynamics of the network; e) Trajectory 
starting from A?, (at t = 0 the total mass is in A3), undergoing two transitions 
at ti and t 2 - The simulation has been performed for kinetic constants ^ 

with e = 1/50. On top, concentration of species (concentrations of Al,A 4 ,Ag 
are negligible everywhere). At bottom, orders of concentrations (computed as 
log£(a::i)) with continuous lines if species is tropically equilibrated, dotted lines 
if not. 
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3 Tropical equilibrations of nonlinear networks 
with polynomial rate functions 


In this section we consider nonlinear biochemical networks described by mass 
action kinetics 

^ = X! kjSijX°‘\ l<i<n, (9) 

3 

where kj > 0 are kinetic constants, Sij are the entries of the stoichiometric 
matrix (uniformly bounded integers, < s, s is small), aj = (a\,... ,a^) 

are multi-indices, and a;“j = Xi^ .. .Xn", where a{ are positive integers. 

For chemical reaction networks with multiple timescales it is reasonable to 
consider that kinetic parameters have different orders of magnitudes. This can 
be conveniently formalized by considering that parameters of the kinetic models 
([^ can be written as kj = kje'^^ (10). The exponents 7 j are considered to be 
integer or rational. For instance, the approximation 7 ^ = round(log(fcj)/log(e)) 
produces integer exponents, whereas 7 ^- = round(dlog(fcj)/log(e))/c? produces 
rational exponents, where round stands for the closest integer (with half-integers 
rounded to even numbers) and d is a strictly positive integer. Kinetic parameters 
are fixed. In contrast, species orders vary in the concentration space and have 
to be calculated as solutions to the tropical equilibration problem. To this aim, 
the network dynamics is first described by a rescaled ODE system 

3 

where f3j{a) = 7 ^ -I- {a,aj) (12), and (,) stands for the dot product. 

The r.h.s. of each equation in 0 is a sum of multivariate monomials in 
the concentrations. The orders f3j indicate how large are these monomials, in 
absolute value. A monomial of order /ij dominates another monomial of order 

fij' if jij < 11 ji. 

The tropical equilibration problem consists in the equality of the orders of 
at least two monomials one positive and another negative in the differential 
equations of each species. More precisely, we want to find a vector a such that 


min ( 7 , 

hSij>0 


+ (a, a 


,)) = 


ST3 


+ (a, a 


;)) 


(13) 


j,Sij<0' 

Computing tropical equilibrations from the orders of magnitude of the model 
parameters is a NP-complete problem, cf. m- However, methods based on the 
Newton polytope |15j or constraint logic programming |16j exploit the sparseness 
and redundance of the system to effectively obtain sets of solutions. The equa¬ 
tion! 13) is related to the notion of tropical hypersurface. A tropical hypersurface 
is the set of vectors a G M” such that the minimun min^- 5 .^.^o( 7 j + (o-jCtj)) is 
attained for at least two different indices j (with no sign conditions). Tropi¬ 
cal prevarieties are finite intersections of tropical hypersurfaces. Therefore, our 
tropical equilibrations are subsets of tropical preverieties. The sign condition 
in (13) was imposed because species concentrations are real positive numbers. 
Compensation of a sum of positive monomials is not possible for real values of 
the variables. 

Species timescales. The timescale of a variable xt is given by A ^ = A ^ 
whose order is Vi = min{/rj|S'ij ^ 0} — (14). The order vi indicates how fast 

is the variable xi (if Vii < Ui then Xi/ is faster than Xi) . 


6 




Partial tropical equilibrations. It is useful to extend the tropical equilibration 
problem to partial equilibrations, that means solving (131 only for a subset 
of species. This is justified by the fact that slow species do not need to be 
equilibrated. In order to have a self-consistent calculation we compute the 
species timescales by (141. A partial equilibration is consistent < v for all 
non-equilibrated species i. > 0 is an arbitrarily chosen threshold indicating 
the timescale of interest. 

Tropical equilibrations, slow invariant manifolds and metastable states. In 
dissipative systems, fast variables relax rapidly to some low dimensional attrac¬ 
tive manifold called invariant manifold [3] that carries the slow mode dynamics. 
A projection of dynamical equations onto this manifold provides the reduced 
dynamics [5]. This simple picture can be complexified to cope with hierarchies 
of invariant manifolds and with phenomena such as transverse instability, ex¬ 
citability and itineracy. Firstly, the relaxation towards an attractor can have 
several stages, each with its own invariant manifold. During relaxation towards 
the attractor, invariant manifolds are usually embedded one into another (there 
is a decrease of dimensionality) [2]. Secondly, invariant manifolds can lose local 
stability, which allow the trajectories to perform large phase space excursions 
before returning in a different place on the same invariant manifold or on a 
different one [7]. We showed elsewhere that tropical equilibrations can be used 
to approximate invariant manifolds for systems of polynomial differential equa¬ 
tions [siiiniiii]- Indeed, tropical equilibration are defined by the equality of 
dominant forces acting on the system. The remaining weak non-compensated 
forces ensure the slow dynamics on the invariant manifold. Tropical equilibra¬ 
tions are thus different from steady states, in that there is a slow dynamics. In 
this paper we will use them as proxies for metastable states. 

Branches of tropical equilibrations and connectivity graph. For each equation 
i, let us define Mi{a) = argmin(p,j(a), > 0) = argmin(/ij (a), S'ij < 0) (15), 

j j 

in other words Mi denotes the set of monomials having the same minimal order 
pLi. We call tropically truncated system the system obtained by pruning the 
system (111, i.e. by keeping only the dominating monomials. 

( 16 ) 

j&Miia) 


dt 


The tropical truncated system is uniquely determined by the index sets Mi{a), 
therefore by the tropical equilibration a. Reciprocally, two tropical equilibra¬ 
tions can have the same index sets Mi{a) and truncated systems. We say that 
two tropical equilibrations oi, 02 are equivalent iff Mi{ai) = Mi(a 2 ),for all i. 
Equivalence classes of tropical equilibrations are called branches. A branch B 
with an index set Mi is minimal if M' C Mi for all i where M' is the index 
set B' implies B' = B or B' = 0. Closures of equilibration branches are de- 
hned by a hnite set of linear inequalities, which means that they are polyhedral 
complexes. Minimal branches correspond to maximal dimension faces of the 
polyhedral complex. The incidence relations between the maximal dimension 
faces (n — I dimensional faces, where n is the number of variables) of the poly¬ 
hedral complex define the conneetivity graph. More precisely, minimal branches 
are the vertices of this graph. Two minimal branches are connected if the cor¬ 
responding faces of the polyhedral complex share a n — 2 dimensional face. In 
terms of index sets, two minimal branches with index sets M and M' are con- 
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nected if there is an index set M" such that M' C M” and Mi C M" for all 
i. 


Tropical equilibrations and monomolecular networks. Eqs.( 13) have a simpler 
form in the case of monomolecular networks 


min (7y+ai)= min ( 7 ^-,+ a*) (17) 

_7GPred(2) jGoucc(2) 

where Pred(*) = {j\{j,i) £ A}, Succ(i) = {j\{i,j) £ A} are the sets of prede¬ 
cessors and successors of the node i in the digraph G. 

Let us recall that by min-plus algebra we understand the semi-ring (R U 
{oo}, 0 ,(g)) where the two operations are defined as x (B y = min{a;,j/} and 
X ® y = X + y. In other words the addition and the min operation play the role 
of min-plus multiplication and addition, respectively. Therefore Eqs. 0 are 
linear in the unknowns ai. Computing tropical equilibrations of monomolecular 
networks boils down to solving linear equations in min-plus algebra. For linear 
tropical systems there are fast algorithms [HIS]. 

We have tested the tropical equilibration conditions 0 for the trajectories 
of the monomolecular network presented in Figure by checking if the abso¬ 
lute value of the difference between the r.h.s and l.h.s of 0 is smaller than 
a threshold. The result is illustrated in Fig. 0 - For this model, the tropical 
equilibration solutions are changing along the trajectory. This can been seen by 
following the orders of the concentrations along the trajectories. These orders 
change by integers at transition points. Furthermore, at transition points some 
of the variables that where not previously equilibrated, become equilibrated. 
The analysis of the tropical equilibrations finds the transitions previously de¬ 
tected in Section]^ from the approximated eigenvalues and eigenvectors (ti and 
t 2 for this example) but adds some more. For instance, species A1 equilibrates 
at the timescale 1 /ki = 10. This was not taken into account in the descrip¬ 
tion of the automaton in Figure Id) because the species A1 is fast and can not 
accumulate. 


4 Learning a finite state machine from a nonlin¬ 
ear biochemical network 

We are using the algorithm based on constraint solving introduced in |16j to 
obtain all rational tropical equilibration solutions a = (oi, 02 ,..., a„) within 
a box \ai\ < b, b > 0 and with denominators smaller than a hxed value d, 
0 = Pi/q, Pi,q are positive integers, q < d. The output of the algorithm is 
a matrix containing all the tropical equilibrations within the defined bounds. 
A post-processing treatment is applied to this output consisting in computing 
truncated systems, index sets, and minimal branches. Tropical equilibrations 
minimal branches are stored as matrices Ai, A 2 ,..., A;,, whose lines are tropical 
solutions within the same branch. Here b is the number of minimal branches. 

Our method computes numerical approximations of the tropical prevariety. 
Given a value of e, this approximation is better when the denominator bound d is 
high. At fixed d, the dependence of the precision on e follows more intricate rules 
dictated by Diophantine approximations. For this reason, we systematically test 
that the number b and the truncated systems corresponding to minimal branches 
are robust when changing the value of e. 



Trajectories x{t) = {xi{t),..., Xn{t)) of the smooth dynamical system are 
generated with different initial conditions, chosen uniformly and satisfying the 
conservation laws, if any. For each time t, we compute the Euclidian distance 
di{t) = min^g^^ Wv ~ ^oge{x{t))\\ , where ||*|| denotes the Euclidean norm and 
logg(a;) = (logxi/log(e),...,logx„/log(e)). This distance classifies all points 
of the trajectory as belonging to a tropical minimal branch. The result is a 
symbolic trajectory Si, S 2 , ■ ■ ■ where the symbols Si belong to the set of minimal 
branches. In order to include the possibility of transition regions we include 
an unique symbol t to represent the situations when the minimal distance is 
larger than a fixed threshold. We also store the residence times ti,T 2 , ... that 
represent the time spent in each of the state. 

The stochastic automaton is learned as a homogenous, finite states, contin¬ 
uous time Markov process, defined by the lifetime (mean sojourn time) of each 
state Ti, 1 < i < b and by the transition probabilities pij from a state i to 
another state j. We use the following estimators for the lifetimes and for the 
transition probabilities: 

Ti = (18) 

n n 

P'i-d ~ ll-Sn=i.Sn+l=j)/(^^ ll-S„=i)) * J (l^) 

n n 

As a case study we consider a nonlinear model of dynamic regulation of Trans¬ 
forming Growth Eactor beta TGF-/3 signaling pathway proposed in [T]. This 
model has a dynamics defined by n = 18 polynomial differential equations and 
25 biochemical reactions. The paper [T] proposes three versions of the mecha¬ 
nism of interaction of TIFI 7 (Transcriptional Intermediary Factor 1 7 ) with the 
Smad-dependent TGF-/3 signaling. We consider here the version in which TIFl 
interacts with the phosphorylated Smad2-Smad4 complexes leading to dissoci¬ 
ation of the complex and degradation of Smad4. The results are similar for the 
other versions of this model. The example was chosen because it is a medium 
size model based on polynomial differential equations. The computation of the 
tropical equilibrations for this model shows that there are 9 minimal branches 
of full equilibrations (in these tropical solutions all variables are equilibrated). 
The connectivity graph of these branches and the learned automaton are shown 
in Figure]^ The study of this example shows that branches of tropical equilibra¬ 
tion can change on trajectories of the dynamical system. Furthermore, all the 
observed transitions between branches are contained in the connectivity graph 
resulting from the polyhedral complex of the tropical equilibration branches. 

The transition probabilities of the automaton are coarse grained properties 
of the statistical ensemble of trajectories for different initial conditions. Given a 
state and a minimal branch close to it, it will depend on the actual trajectory to 
which other branch the system will be close to next. However, when initial data 
and the full trajectory are not known, the automaton will provide estimates of 
where we go next and with which probability. For the example studied, the 
branch B1 is a globally attractive sink: starting from anywhere, the automa¬ 
ton will reach B1 with probability one. This branch contains the unique stable 
steady state of the initial model. Figure [^bottom right shows the structure of 
most probable branches, the ones in which the systems spends most of his time. 
The branches Bl, B3 and B2 correspond to different compositions of the mem¬ 
brane and of the endosome, rich in the receptor RI, rich in the receptor RII and 
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Figure 2: TGF/3 model. Upper left: Connectivity graph of tropical minimal 
branches; upper right: finite state automaton; bottom left: trajectories with 
jumps and distances to minimal branches; the closest branch changes with time 
along the trajectory; bottom right: first three tropical equilibrations minimal 
branches in various projections in concentration orders space. The variables 
RI, RII, LR are membrane receptors (signaling input layers) concentration or¬ 
ders, whereas pS2n, S4n, pS24n are nuclear transcription factors and complexes 
(effectors) concentration orders. The structure tropical branches shows that 
composition of input layers is more flexible (varies on planes) than the concen¬ 
trations of effectors (vary on lines). 


rich in both types of receptors, respectively. Even if this compositon is changed 
on wide domains of orders (planes in the space of orders), the concentrations of 
effectors are robust (are more constained than the concentrations of receptors). 

5 Conclusion 

We have presented a method to coarse grain the dynamics of a smooth biochemi¬ 
cal reaction network to a discrete symbolic dynamics of a finite state automaton. 
The coarse graining was obtained by two methods, approximated eigenvectors 
for mono-molecular networks and minimal branches of tropical equilibrations 
for more general mass action nonlinear networks. The two methods are com¬ 
patible one to another, because when applied to monomolecular networks the 
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method based on tropical geometry detects all the transitions indicated by ap¬ 
proximated eigenvectors. For both methods the automaton has a small number 
of states, less than the number of species in the first method and the number of 
minimal tropical branches in the second method. The coarse grained automa¬ 
ton can be used for studying statistical properties of biochemical networks such 
as occurrence and stability of temporal patterns, recurrence, periodicity and 
attainability problems. The coarse graining can be performed in a hierarchical 
way. For the nonlinear example studied in the paper we computed only the full 
tropical equilibrations that stand for the lowest order in the hierarchy (coarsest 
model). As discussed in Section 3 we can also consider partial equilibrations 
when slow variables are not equilibrated and thus refine the automaton. Our 
approach extends the notion of steady states of a network and propose a sim¬ 
ple recipe to characterize and detect metastable states. Most likely metastable 
states have biological importance because the network spends most of its time 
in these states. The itinerancy of the network, described as the possibility of 
transitions from one metastable state to another is paramount to the way neural 
networks compute, retrieve and use information |18j and can have similar role 
in biochemical networks. 
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Appendix 1: Proof of Proposition 


2.2 


Let us consider that = 1. 


Taking = 0 for all predecessors j of k and for all other nodes that lead to k 
by the flow d) satisfy Eq.(|^(main body text) with A = —Kk- The same is valid 
for all the nodes that do not lead to k and are not accessible from k. Remain 
the nodes that are accessible from k. Let j be such a node. Then j = 
for some m > 0. Eq.(|^(main body text) implies that 

= (-fik + for 1 < Z < m. 


Thus r| 






X ... X 


hfc) 




. Suppose that Kk < 


(fc) -Kk+Kil,(k) " -Kfc+K'.j2(fc) 

K^Hk) for Z = 1,..., TO — 1 and K^mf^k) < Using Lemm? 2.1 main body text) 
= — 1. If any of the previous inequality does not hold then at 


it follows 


qfe) 


least one factor in the expression of vanishes and the remaining factors 

are finite, thus = 0. Consider now that Zf = 1. Taking li = 0 for all the 

nodes j that can be obtained from k and for all other nodes that do not lead to 
k by the flow $ satisfy Eq.Q(main body text) with A = —Kk- The remaining 
nodes are all leading to k. Let j be such a node. Then k = $'"(j) for some 
TO > 0. Eq.Q (main body text) implies that 

K$<- 1 (j)Z|iq) = {-Kk + K$i-i(j))^$i-i(j), for 1 < Z < TO. 


TTpncp / • = -^- 

nence 


for all Z = 0, 






X ... X 


pj) 






Suppose that «$!(_,) > K-k, 

, TO—1. Using Lemms 2.1 main body text) it follows Ij = 1. If one 

, TO — 1 then the corresponding 


of these inequalities is not satisfied for a Z = 0, 
factor in the expression of vanishes and 1^=0. 

The above formulas cover the zero eigenvalue case if we consider that Kk = 0 
for k being the sink. It follows that = 1 and rj=0 elsewhere. Furthermore, 
Z° = 1 for all j. 


Appendix 2: Algorithm for reduction of monomolecular networks 
with total separation separation. This algorithm consists of three steps. 

I. Constructing of an auxiliary reaction network: pruning. 

For each Ai branching node (substrate of several reactions) let us define Ki 
as the maximal kinetic constant for reactions Ai —>■ Aj: Ki = maxjjfcji}. For 
correspondent j we use the notation j = ^i(f): (j){i) = a.i'gmaxj{kji}. 

An auxiliary reaction network V is the set of reactions obtained by keeping 
only Ai —^ with kinetic constants Ki and discarding the other, slower 

reactions. Auxiliary networks have no branching, but they can have cycles and 
confluences. The correspondent kinetic equation is 

C-i — KiCi T ^ ^ 

If the auxiliary network contains no cycles, the algorithm stops here. 

II gluing cycles and restoring cycle exit reactions 

In general, the auxiliary network V has several cycles Ci, (72 ,... with periods 
Ti,T2 , ... > I. 

These cycles will be “glued” into points and all nodes in the cycle Ci, will 
be replaced by a single vertex A*. Also, some of the reactions that were pruned 
in the first part of the algorithm are restored with renormalized rate constants. 
Indeed, reaction exiting a cycle are needed to render the correct dynamics: 
without them, the total mass of the cycle is conserved, with them the mass 


I 











can also slowly leave the cycle. Reactions A ^ B exiting from cycles (A € Ci, 
B ^ Ci) are changed into ^ B with the rate constant renormalization: let 
the cycle C* be the following sequence of reactions —)■ ^2 —> ■■■Ar- —)■ Ai, 
and the reaction rate constant for Ai —)■ is ki {kn for At- —>■ ^ 1 ). For the 
limiting (slowest) reaction of the cycle Ci we use notation kn^n i- li A = Aj and 
k is the rate reaction for A ^ B, then the new reaction A^ ^ B has the rate 
constant kkumi/kj. This rate is obtained using quasi-stationary distribution 
for the cycle. If kinetic constants are expressed as powers of a small positive 
parameter e, i.e., ii k = C, then the order of the constant has to be changed 
according to the rule 7 —>■ ') + ^um. — Iji where 7 , 'yu-mi, Ij are the orders of the 
constants k, k^jn i and kj, respectively. 

The new auxiliary network is computed for the network of glued cycles. 
Then we decompose it into cycles, glue them, iterate until a acyclic network is 
obtained V”. 

Ill Restoring cycles 

The dynamics of species inside glued cycles is lost after the second part. A 
full multi-scale approximation (including relaxation inside cycles) can be ob¬ 
tained by restoration of cycles. This is done starting from the acyclic auxiliary 
network V" back to through the hierarchy of cycles. Each cycle is restored 
according to the following procedure: 

For each glued cycle node A™, node of V™, 

• Recall its nodes —)• A^~^ —>• —)• they form a cycle of 

length Ti- 

• Let us assume that the limiting step in A™ is A™“^ —)• A”(“^ 

• Remove A™ from V"* 

• Add n vertices A™“\...A™“^ to V™ 

• Add to V™ reactions A”(“^ —>• A)^~^ —>• ...A™“^ (that are the cycle reac¬ 
tions without the limiting step) with correspondent constants from 

• If there exists an outgoing reaction A™ B in V™ then we substitute it 
by the reaction A™“^ —> B with the same constant, i.e. outgoing reactions 
A™ —> ... are reattached to the beginning of the limiting steps 

• If there exists an incoming reaction in the form B —)■ A™, hnd its prototype 

in and restore it in V™ 

• If in the initial V'" there existed a “between-cycles” reaction A™ —)■ A™ 

then we find the prototype in A —)• R, and substitute the reaction 

by A™”^ —)• B with the same constant, as for A™ —>• A™ (again, the 
beginning of the arrow is reattached to the head of the limiting step in 

AT) 
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Figure 3: The successive steps of the reduction algorithm, illustrated for the 
prism model used in the paper, a) is the initial model; b) is the auxiliary net¬ 
work resulting from step I, pruning; c) is the result of gluing 3 species cycles 
and renormalizing the exit reactions (the constants of orders 3, 7,10,8 are renor¬ 
malized to 3 -I- 6 — 1 = 8,7 -I- 6 — 6 = 7, 10 -I- 6 — 4 = 12, and 8 4-9 — 2 = 15, 
respectively); d) is the auxiliary network after one more iteration; e) results 
from gluing and then restoring the 3 species cycles without the limiting step 
(constant of order 15); f) results from restoring the single species cycles without 
their limiting steps. 
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Appendix 3: Description of the TGFb model used in this paper. The 

model is described by the following system of differential equations 
dxi 


dt 

dX2 

dt 

dX3 

dt 

da;4 

dt 

dxs 

dt 

dxe 

dt 

dx'j 

dt 

dxs 

dt 

dxg 

dt 

dxio 

dt 

dxii 

dt 

da:i 2 

dt 

dxia 

dt 

da; 14 
dt 

da:i5 

dt 

da:i6 

dt 

da:i7 

dt 

da:i 8 

dt 


= k2X2 — kixi — kiQXiXii 
= kixx — k2X2 + knk34Xe 
= k^Xi - ksxs + krxr + ks^ksrxis - k^x^x^ 

= ksxs - kzXA + kgxs - ksX4Xe 

= k^XQ — k4X5 + k^XY + 2kiiXQ — 2kioxl — kgx^x^ + kigXiXn 
= k4X5 - k^XQ + kgXs + 2fci3a:io - 2ki2xl - knkg4XQ + ksikgtjXs - ksX4Xe 
= keX3X5 - X7{kY + ku) 

= ki4XY - kgXs - ksikggXs + ksX4Xg 
= kigx\ — Xgikii + tcis) 

= ^152^9 ~ ^132;i0 + ki2x\ 

= ^232^14 ~ fc302;il 

= ^18 ~ a;i2(fc20 + ^26) + fc302;il + fc272;i5 ~ kggks^XigXig 

= kig — a;i3(fc21 + k2s) + fc302;il + fc29a;i6 ~ k 22 kzY,Xi 2 Xig 
= k22kggXi2Xig — Xi4(fc23 + ^24 + ^ 25 ) 

= k2gXi2 — kgYXiY, 

= k2sXi3 — kggXiQ 

= kgikgQXs — kg2Xn 
= kggXn — kggk^YXis 


These variables are as follows: 

• Receptors on membrane: a;i 2 = RI, a;i 3 = RII, a;i 4 = LR. 

• Receptors in the endosome: a;ii = LRe, X 15 = RIe, Xiq = Rile. 

• Transcription factors and complexes in cytosol: a;i = S2c, 0:3 = S4c, 0:5 = 
pS2c, xy = pS24c, xg = pS22c, xis = S4ubc. 

• Transcription factors and complexes in the nucleus: xg = S2n, 0:4 = S4n, 
a ;6 = pS2n, a ;8 = pS24n, a;io = pS22n, a;i 7 = S4ubn. 
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